Scenario

Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.

Load packages

library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement 
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed 
library(bayesplot) # needed for MCMC diagnostics 
library(DHARMa) # model validation 
library(ggdist) # partial plots 
library(tidybayes) # partial plots 
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans 
library(ggpubr) # arranging figures 
library(mgcv) # GAM 
library(modelr) # plotting
library(segmented) # piece wise regression 
library(strucchange) #piece wise regression

Set working directory

knitr::opts_knit$set(root.dir=working.dir)

Import data

This data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty

lat_resp_dat <- read_delim("C:/Users/jc527762/OneDrive - James Cook University/PhD dissertation/Data/Chapter4_Latitude/import_files/lat_resp_dat_no_outliers.txt", 
    delim = "\t", escape_double = FALSE, 
    col_types = cols(area = col_skip(), rate.a.spec = col_skip()), 
    trim_ws = TRUE)

Data manipulation

rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |> 
  mutate(dev.temp = as.factor(dev.temp), 
         replicate = str_sub(sampleID, -1,-1), 
         population = factor(population)) |>
                                           # number of observations = 5758
  filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
         sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
         sampleID != "60.LCKM.152.30.1"  # 5637 - 76 = 5542
  ) |> 
  filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes) 
  group_by(sampleID) |> 
  mutate(max_value_index = which.max(rate.output), 
         row_number = row_number(), 
         max.rate = max(rate.output), 
         low.rate = mean(rate.output[order(rate.output)[1:10]]), 
         net.rate = max.rate - low.rate) |> 
  filter(row_number <= max_value_index) |>
  ungroup() |> 
  mutate(region = factor(region), 
         dev.temp =factor(dev.temp), 
         level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
                           tank >= 200 & tank <= 299 ~ 2,
                           tank >= 300 & tank <= 399 ~ 3,
                           TRUE ~ NA_real_)), 
         female = factor(female)) |> 
  unite("dev.temp.region",c(dev.temp,region), remove=FALSE) |> 
  mutate(dev.temp.region = factor(dev.temp.region))

Exploratory data analysis

2nd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 2, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

3rd order polynomial

eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) + 
  geom_point(alpha=0.5, color = "black") + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  theme_classic() +
  ggtitle("All data points - 2nd order polynomial")

eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ region) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial") 

eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ dev.temp) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Region - 2nd order polynomial")

eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) + 
  geom_point(alpha=0.5) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE), color = "red") + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) + 
  geom_point(alpha=0.1) + 
  geom_smooth(method="lm", 
              se=TRUE, 
              fill=NA,
              formula = y~poly(x, 3, raw=TRUE)) + 
  facet_wrap(~ population) + theme_classic() +
  theme(legend.position = "none") +
  ggtitle("Population - 2nd order polynomial")

eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
          nrow = 6, 
          ncol=1); eda.fig

Fit model [random factors]

First we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.

Hypothesis test will include:

  1. Null model
  2. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK)
  3. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL)
  4. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)
  5. EGG_COUNT ~ 1 + (1| FEMALE) + (1| TANK) + (1| LEVEL) + (1| POPULATION)+ (1| CLUTCH_ORDER)
lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output))) 
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

rate.priors2 <- prior(normal(0.635, 0.20), class="Intercept") +
  prior(student_t(3, 0, 0.45), class = "sigma")
rate.priors <- prior(normal(430, 0.25), class="Intercept") +  
  prior(student_t(3, 0, 195), class = "sigma")

Models

Null

f.model.null <- bf(rate.output2 ~ 1, 
                   family = gaussian()) 

model.null <- brm(f.model.null,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model.null, file="./03.CTmax_resp_data_analsyis_files/model.null.RDS")

Model1

f.model1 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank), 
                          family=gaussian())

model1 <- brm(f.model1,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1, file="./03.CTmax_resp_data_analsyis_files/model1.RDS")

Model2

f.model2 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank) + (1| level), 
                          family=gaussian())

model2 <- brm(f.model2,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 26 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 753 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model2, file="./03.CTmax_resp_data_analsyis_files/model2.RDS")

Model3

f.model3 <- bf(rate.output2 ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population), 
               family=gaussian())

model3 <- brm(f.model3,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 18 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 85 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model3, file="./03.CTmax_resp_data_analsyis_files/model3.RDS")

Model4

f.model4 <- bf(rate.output2 ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order), 
                        family=gaussian())

model4 <- brm(f.model4,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE),  
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 74 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 14 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model4, file="./03.CTmax_resp_data_analsyis_files/model4.RDS")

LOO

LOO(model.null,model1, model2,model3,model4) 
## Output of model 'model.null':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -1446.4  81.9
## p_loo         3.9   0.3
## looic      2892.8 163.8
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [1.0, 1.0]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo     91.4  83.1
## p_loo        56.1   2.3
## looic      -182.8 166.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo     91.3  83.2
## p_loo        56.3   2.3
## looic      -182.7 166.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo     91.3  83.1
## p_loo        56.3   2.3
## looic      -182.5 166.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model4':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo     90.8  83.2
## p_loo        56.8   2.3
## looic      -181.7 166.4
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##            elpd_diff se_diff
## model1         0.0       0.0
## model2        -0.1       0.4
## model3        -0.1       0.4
## model4        -0.5       0.5
## model.null -1537.8      77.4

Fit model [fixed factors - polynomials]

Prior investigation

lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output2)), 
            sd(na.omit(rate.output2))) ; lat_resp_dat2 |> 
  group_by(region,dev.temp) |> 
  summarise(mean(na.omit(rate.output)), 
            sd(na.omit(rate.output))) 
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

Priors

# change to reflect rate.output2
rate.priors2 <- prior(normal(0.635, 0.20), class="Intercept") + 
  prior(normal(0, 0.20), class="b") +
  prior(student_t(3, 0, 0.45), class = "sigma")
  
  rate.priors <- prior(normal(430, 0.25), class="Intercept") + 
  prior(normal(0, 40), class="b") +
  prior(student_t(3, 0, 195), class = "sigma")

Models

Model1.p1

Model1.p1

f.model1.p1 <- bf(rate.output2 ~ 1 + 
                         dev.temp*region + Temperature +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p1 <- brm(f.model1.p1,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p1, file="./03.CTmax_resp_data_analsyis_files/model1.p1.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 972987, p-value < 2.2e-16
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p1$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p1$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p1$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p1$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p1$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

model1.p2

model1.p2

f.model1.p2 <- bf(rate.output2 ~ 1 + 
                         dev.temp*region + poly(Temperature, 2) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p2 <- brm(f.model1.p2,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p2, file="./03.CTmax_resp_data_analsyis_files/model1.p2.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1044481, p-value < 2.2e-16
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p2$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p2$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p2$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p2$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p2$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

model1.p3

model1.p3

f.model1.p3 <- bf(rate.output2 ~ 1 + 
                         dev.temp*region + poly(Temperature, 3) +
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank) , 
                         family=gaussian()) 

model1.p3 <- brm(f.model1.p3,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(model1.p3, file="./03.CTmax_resp_data_analsyis_files/model1.p3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1078090, p-value < 2.2e-16
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.p3$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.p3$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.p3$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.p3$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.p3$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

LOO

LOO(model1.p1, model1.p2,model1.p3) 
## Output of model 'model1.p1':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo   1452.7  84.8
## p_loo        60.6   2.7
## looic     -2905.4 169.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p2':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo   1628.3  88.7
## p_loo        60.9   2.8
## looic     -3256.6 177.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model1.p3':
## 
## Computed from 1800 by 4776 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo   1701.8  84.3
## p_loo        61.8   2.7
## looic     -3403.7 168.7
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
## 
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##           elpd_diff se_diff
## model1.p3    0.0       0.0 
## model1.p2  -73.5       9.2 
## model1.p1 -249.1      21.5

Fit model [fixed factors - GAM]

GAM model k=4

GAM model k=4

f.model1.gamk4 <- bf(rate.output2 ~ 1 + 
                         t2(dev.temp,region,Temperature, k=4, bs="fs", by=dev.temp.region) + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk4 <- brm(f.model1.gamk4,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 2 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 899 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model1.gamk4, file="./03.CTmax_resp_data_analsyis_files/model1.gamk4.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="dev.temp.region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1199428, p-value < 2.2e-16
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk4$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk4$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk4$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk4$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk4$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

### {-}

GAM model k=3

GAM model k=3

f.model1.gamk3 <- bf(rate.output2 ~ 1 + 
                         t2(dev.temp,region,Temperature, k=3, bs="fs", by=dev.temp.region) + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

model1.gamk3 <- brm(f.model1.gamk3,
              data = lat_resp_dat2, 
              prior = rate.priors2,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 25 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 1422 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(model1.gamk3, file="./03.CTmax_resp_data_analsyis_files/model1.gamk3.RDS")

DHARMa

#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")

#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds), 
                              observedResponse = lat_resp_dat2$rate.output, 
                              fittedPredictedResponse = apply(preds, 2, mean), 
                              integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1160466, p-value < 2.2e-16
## alternative hypothesis: two.sided

MCMC diagnostics

stan_trace(model1.gamk3$fit)
## 'pars' not specified. Showing first 10 parameters by default.

stan_ac(model1.gamk3$fit) 
## 'pars' not specified. Showing first 10 parameters by default.

stan_rhat(model1.gamk3$fit) 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_ess(model1.gamk3$fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

stan_dens(model1.gamk3$fit, separate_chains = TRUE)
## 'pars' not specified. Showing first 10 parameters by default.

### {-}

AIC(model1.gamk3, model1.gamk4)
## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observations

Partial effects plots

conditional_effects

model1.gamk4 |> 
  conditional_effects(spaghetti = TRUE, ndraws=200) 

Model investigation

summary

summary(model1.gamk4)
## Warning: There were 2 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rate.output2 ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs", by = dev.temp.region) + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank) 
##    Data: lat_resp_dat2 (Number of observations: 4776) 
##   Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
##          total post-warmup draws = 1800
## 
## Smooth Terms: 
##                                                             Estimate Est.Error
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)        1.49      1.26
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)        1.52      1.26
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)     2.33      1.37
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)     2.27      1.33
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)        1.94      1.12
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)        1.99      1.25
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)     2.17      1.25
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)     2.26      1.27
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)        1.42      1.06
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)        1.54      1.23
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)     1.70      1.00
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)     1.71      1.21
##                                                             l-95% CI u-95% CI
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)        0.46     4.10
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)        0.43     4.44
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)     0.86     5.70
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)     0.88     5.63
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)        0.70     4.77
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)        0.75     5.15
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)     0.79     5.35
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)     0.85     5.61
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)        0.41     4.36
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)        0.46     4.50
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)     0.58     4.19
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)     0.55     4.76
##                                                             Rhat Bulk_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)    1.00     1713
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)    1.00     1567
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1.00     1550
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1.00     1819
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)    1.00     1616
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)    1.00     1855
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1.00     1561
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1.00     1620
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)    1.00     1586
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)    1.00     1787
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1.00     1663
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1.00     1780
##                                                             Tail_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1)        1631
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2)        1503
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1)     1630
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2)     1872
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1)        1536
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2)        1709
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1)     1610
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2)     1687
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1)        1792
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2)        1743
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1)     1477
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2)     1591
## 
## Group-Level Effects: 
## ~female (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.10      0.04     0.03     0.20 1.00      878     1032
## 
## ~tank (Number of levels: 52) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.02     0.11     0.18 1.00     1499     1454
## 
## Population-Level Effects: 
##                                  Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                            0.49      0.08     0.32     0.64 1.00
## scalemasscenterEQTRUEscaleEQTRUE     0.26      0.01     0.25     0.27 1.00
##                                  Bulk_ESS Tail_ESS
## Intercept                            1635     1644
## scalemasscenterEQTRUEscaleEQTRUE     1750     1717
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.16      0.00     0.16     0.16 1.00     1883     1743
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#car::Anova(model1.geo.gamma)

R2

model1.gamk4 |> bayes_R2(summary = FALSE) |> median_hdci()

tidyMCMC

tidyMCMC(model1.gamk4, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

gather draws

#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |> 
  median_hdci()

bayesplot

model1.gamk4 |> mcmc_plot(type='intervals')

emmeans - pariwise

model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)
## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()
## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> regrid() 
## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)
##  region  dev.temp emmean lower.HPD upper.HPD
##  core    28        0.531     0.412     0.659
##  leading 28        0.550     0.422     0.698
##  core    30        0.635     0.526     0.737
##  leading 30        0.498     0.352     0.640
##  core    31        0.618     0.494     0.725
##  leading 31        0.607     0.487     0.730
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

emtrends - pariwise

emtrends(model1.gamk4, specs=pairwise~dev.temp*region, var = "Temperature")  
## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)
## $emtrends
##  dev.temp region  Temperature.trend lower.HPD upper.HPD
##  28       core              0.05589    0.0364   0.07286
##  30       core              0.00379   -0.0142   0.01943
##  31       core              0.02995    0.0132   0.04545
##  28       leading          -0.01301   -0.0325   0.00544
##  30       leading           0.00435   -0.0153   0.02521
##  31       leading           0.01009   -0.0071   0.02705
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast                                 estimate lower.HPD upper.HPD
##  dev.temp28 core - dev.temp30 core        0.052097  0.025961   0.07356
##  dev.temp28 core - dev.temp31 core        0.025937  0.000936   0.04901
##  dev.temp28 core - dev.temp28 leading     0.068667  0.043783   0.09549
##  dev.temp28 core - dev.temp30 leading     0.051514  0.025472   0.07994
##  dev.temp28 core - dev.temp31 leading     0.046443  0.020681   0.06957
##  dev.temp30 core - dev.temp31 core       -0.026011 -0.049167  -0.00357
##  dev.temp30 core - dev.temp28 leading     0.017010 -0.007222   0.04247
##  dev.temp30 core - dev.temp30 leading    -0.000309 -0.025766   0.02605
##  dev.temp30 core - dev.temp31 leading    -0.005979 -0.031522   0.01642
##  dev.temp31 core - dev.temp28 leading     0.042749  0.019082   0.06782
##  dev.temp31 core - dev.temp30 leading     0.025722 -0.000808   0.05269
##  dev.temp31 core - dev.temp31 leading     0.019962 -0.004044   0.04130
##  dev.temp28 leading - dev.temp30 leading -0.017101 -0.046804   0.00952
##  dev.temp28 leading - dev.temp31 leading -0.023032 -0.049326   0.00121
##  dev.temp30 leading - dev.temp31 leading -0.005859 -0.032276   0.02083
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

probabilities

mtsqst <- model1.gamk4 |> emmeans(pairwise ~ region*dev.temp)
## NOTE: A nesting structure was detected in the fitted model:
##     dev.temp.region %in% (dev.temp*region)
mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())

values

values.df <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |> 
  filter(Temperature < 34) |> 
  group_by(dev.temp.region) |> 
  summarise(oxygen_uptake_values = median(.value)) |>
  ungroup() |> 
  separate(dev.temp.region, c("temperature","region"), remove=FALSE); values.df
## Warning: `fitted_draws` and `add_fitted_draws` are deprecated as their names were confusing.
## - Use [add_]epred_draws() to get the expectation of the posterior predictive.
## - Use [add_]linpred_draws() to get the distribution of the linear predictor.
## - For example, you used [add_]fitted_draws(..., scale = "response"), which
##   means you most likely want [add_]epred_draws(...).
## NOTE: When updating to the new functions, note that the `model` parameter is now
##   named `object` and the `n` parameter is now named `ndraws`.

summary figure

resp.plot <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~dev.temp); resp.plot

resp.plot2 <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic() + 
  facet_wrap(~region); resp.plot2

resp.plot3 <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic(); resp.plot3

Pre-inflection point

create data frame

Split posterior of draws from main model into data prior to 33.5 C (earliest inflection point)

pre_infl <- lat_resp_dat2 |> 
  filter(t < 33.5) 

priors

pre_infl |> 
  group_by(dev.temp.region) |> 
  summarise(mean(na.omit(rate.output2)), 
            sd(na.omit(rate.output2)))
pre.priors <- prior(normal(0.57, 0.2), class="Intercept") + 
  prior(normal(0, 0.25), class="b") +
  prior(student_t(3, 0, 0.40), class = "sigma")

linear model

f.pre.model1 <- bf(rate.output2 ~ 1 + 
                         dev.temp*region*Temperature + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

pre.model1 <- brm(f.pre.model1,
              data = pre_infl, 
              prior = pre.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
## Warning: There were 758 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
saveRDS(pre.model1, file="./03.CTmax_resp_data_analsyis_files/pre.model1.RDS")

Model investigation

pre.model1 |> 
  conditional_effects(spaghetti = TRUE, ndraws=200) 

summary

summary(pre.model1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rate.output2 ~ 1 + dev.temp * region * Temperature + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank) 
##    Data: pre_infl (Number of observations: 2439) 
##   Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
##          total post-warmup draws = 1800
## 
## Group-Level Effects: 
## ~female (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.10      0.04     0.03     0.20 1.00     1413     1566
## 
## ~tank (Number of levels: 52) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.02     0.11     0.19 1.00     1410     1336
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                                0.30      0.14     0.03     0.60 1.00
## dev.temp30                              -0.18      0.17    -0.51     0.14 1.00
## dev.temp31                               0.13      0.16    -0.17     0.46 1.00
## regionleading                            0.12      0.16    -0.21     0.42 1.00
## Temperature                              0.01      0.00    -0.00     0.01 1.00
## scalemasscenterEQTRUEscaleEQTRUE         0.25      0.01     0.24     0.26 1.00
## dev.temp30:regionleading                 0.14      0.20    -0.25     0.53 1.00
## dev.temp31:regionleading                 0.14      0.19    -0.23     0.49 1.00
## dev.temp30:Temperature                   0.01      0.01     0.00     0.02 1.00
## dev.temp31:Temperature                  -0.00      0.01    -0.01     0.01 1.00
## regionleading:Temperature               -0.00      0.01    -0.01     0.01 1.00
## dev.temp30:regionleading:Temperature    -0.01      0.01    -0.03     0.00 1.00
## dev.temp31:regionleading:Temperature    -0.01      0.01    -0.02     0.01 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                1993     1737
## dev.temp30                               1959     1738
## dev.temp31                               1846     1679
## regionleading                            1615     1795
## Temperature                              1846     1743
## scalemasscenterEQTRUEscaleEQTRUE         1794     1667
## dev.temp30:regionleading                 1698     1708
## dev.temp31:regionleading                 1776     1635
## dev.temp30:Temperature                   1941     1747
## dev.temp31:Temperature                   1769     1664
## regionleading:Temperature                1565     1648
## dev.temp30:regionleading:Temperature     1723     1683
## dev.temp31:regionleading:Temperature     1738     1678
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.12      0.00     0.12     0.12 1.00     1666     1701
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#car::Anova(model1.geo.gamma)

R2

pre.model1 |> bayes_R2(summary = FALSE) |> median_hdci()

tidyMCMC

tidyMCMC(pre.model1, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

gather draws

#model1.re.wo |> get_variables()
pre.model1 |> gather_draws(`b_.*|sigma`, regex =TRUE) |> 
  median_hdci()

bayesplot

pre.model1 |> mcmc_plot(type='intervals')

emmeans - pariwise

pre.model1 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)
## NOTE: Results may be misleading due to involvement in interactions
pre.model1 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()
## NOTE: Results may be misleading due to involvement in interactions
pre.model1 |> emmeans(pairwise ~ region*dev.temp, type="response") |> regrid() 
## NOTE: Results may be misleading due to involvement in interactions
##  region  dev.temp emmean lower.HPD upper.HPD
##  core    28        0.500     0.370     0.617
##  leading 28        0.563     0.429     0.705
##  core    30        0.650     0.542     0.760
##  leading 30        0.471     0.343     0.626
##  core    31        0.597     0.473     0.708
##  leading 31        0.599     0.462     0.724
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

probabilities

mtsqst <- pre.model1 |> emmeans(pairwise ~ region*dev.temp)
## NOTE: Results may be misleading due to involvement in interactions
mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())

Summary plots

resp.plot <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=33.5, 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(pre.model1, 
                   n=500, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic(); resp.plot

Post-inflection point analysis

post_infl <- lat_resp_dat2 |> 
  filter(t >34.6) 

priors

post_infl |> 
  group_by(dev.temp.region) |> 
  summarise(mean(na.omit(rate.output2)), 
            sd(na.omit(rate.output2)))
post.priors <- prior(normal(0.78, 0.20), class="Intercept") + 
  prior(normal(0, 0.20), class="b") +
  prior(student_t(3, 0, 0.47), class = "sigma")

model

f.post.model1 <- bf(rate.output2 ~ 1 + 
                         dev.temp*region*Temperature + 
                         scale(mass, center=TRUE, scale=TRUE) + 
                         (1| female) + (1| tank), 
                         family=gaussian()) 

post.model1 <- brm(f.post.model1,
              data = post_infl, 
              prior = post.priors,
              warmup = 500, 
              iter = 5000,
              seed=123, 
              cores=2, 
              save_pars = save_pars(all=TRUE), 
              sample_prior = "yes",
              chains = 2, 
              thin = 5, 
              control = list(adapt_delta=0.95)) 
## Compiling Stan program...
## Start sampling
saveRDS(post.model1, file="./03.CTmax_resp_data_analsyis_files/post.model1.RDS")

Model investigation

post.model1 |> 
  conditional_effects(spaghetti = TRUE, ndraws=200) 

summary

summary(post.model1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: rate.output2 ~ 1 + dev.temp * region * Temperature + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank) 
##    Data: post_infl (Number of observations: 1475) 
##   Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
##          total post-warmup draws = 1800
## 
## Group-Level Effects: 
## ~female (Number of levels: 15) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.09      0.05     0.01     0.20 1.00     1111     1422
## 
## ~tank (Number of levels: 52) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.02     0.13     0.22 1.00     1594     1721
## 
## Population-Level Effects: 
##                                      Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                               -7.17      0.32    -7.80    -6.56 1.00
## dev.temp30                              -0.17      0.19    -0.54     0.21 1.00
## dev.temp31                               0.23      0.19    -0.14     0.59 1.00
## regionleading                           -0.18      0.19    -0.56     0.20 1.00
## Temperature                              0.22      0.01     0.21     0.24 1.00
## scalemasscenterEQTRUEscaleEQTRUE         0.28      0.01     0.26     0.30 1.00
## dev.temp30:regionleading                -0.09      0.19    -0.45     0.28 1.00
## dev.temp31:regionleading                 0.02      0.20    -0.36     0.43 1.00
## dev.temp30:Temperature                   0.01      0.01    -0.01     0.02 1.00
## dev.temp31:Temperature                  -0.01      0.01    -0.02     0.01 1.00
## regionleading:Temperature                0.01      0.01    -0.01     0.02 1.00
## dev.temp30:regionleading:Temperature    -0.00      0.01    -0.01     0.01 1.00
## dev.temp31:regionleading:Temperature    -0.00      0.01    -0.01     0.01 1.00
##                                      Bulk_ESS Tail_ESS
## Intercept                                1638     1785
## dev.temp30                               1710     1646
## dev.temp31                               1883     1708
## regionleading                            1735     1697
## Temperature                              1627     1610
## scalemasscenterEQTRUEscaleEQTRUE         1767     1685
## dev.temp30:regionleading                 1624     1708
## dev.temp31:regionleading                 1907     1808
## dev.temp30:Temperature                   1688     1598
## dev.temp31:Temperature                   1901     1613
## regionleading:Temperature                1699     1716
## dev.temp30:regionleading:Temperature     1443     1526
## dev.temp31:regionleading:Temperature     1876     1704
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.20      0.00     0.20     0.21 1.00     1592     1405
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#car::Anova(model1.geo.gamma)

R2

post.model1 |> bayes_R2(summary = FALSE) |> median_hdci()

tidyMCMC

tidyMCMC(post.model1, estimate.method='median', conf.int=TRUE, conf.method='HPDinterval')

gather draws

#model1.re.wo |> get_variables()
post.model1 |> gather_draws(`b_.*|sigma`, regex =TRUE) |> 
  median_hdci()

bayesplot

post.model1 |> mcmc_plot(type='intervals')

emmeans - pariwise

emtrends(post.model1, specs=pairwise~dev.temp*region, var = "Temperature")   
## $emtrends
##  dev.temp region  Temperature.trend lower.HPD upper.HPD
##  28       core                0.222     0.205     0.240
##  30       core                0.228     0.210     0.247
##  31       core                0.217     0.200     0.236
##  28       leading             0.228     0.211     0.246
##  30       leading             0.235     0.212     0.254
##  31       leading             0.222     0.202     0.241
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95 
## 
## $contrasts
##  contrast                                 estimate lower.HPD upper.HPD
##  dev.temp28 core - dev.temp30 core       -0.006634  -0.01665   0.00580
##  dev.temp28 core - dev.temp31 core        0.005451  -0.00620   0.01579
##  dev.temp28 core - dev.temp28 leading    -0.006059  -0.01806   0.00556
##  dev.temp28 core - dev.temp30 leading    -0.012625  -0.03053   0.00557
##  dev.temp28 core - dev.temp31 leading     0.000468  -0.01708   0.01853
##  dev.temp30 core - dev.temp31 core        0.012027  -0.00216   0.02650
##  dev.temp30 core - dev.temp28 leading     0.000288  -0.01667   0.01554
##  dev.temp30 core - dev.temp30 leading    -0.006018  -0.02140   0.00908
##  dev.temp30 core - dev.temp31 leading     0.006744  -0.01352   0.02615
##  dev.temp31 core - dev.temp28 leading    -0.011820  -0.02571   0.00631
##  dev.temp31 core - dev.temp30 leading    -0.017863  -0.04033   0.00123
##  dev.temp31 core - dev.temp31 leading    -0.005220  -0.02090   0.00978
##  dev.temp28 leading - dev.temp30 leading -0.006391  -0.02154   0.00864
##  dev.temp28 leading - dev.temp31 leading  0.006235  -0.00854   0.02151
##  dev.temp30 leading - dev.temp31 leading  0.012649  -0.00721   0.03365
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

probabilities

mtsqst <- post.model1 |> emmeans(pairwise ~ region*dev.temp)
## NOTE: Results may be misleading due to involvement in interactions
mtsqrt2 <- mtsqst$contrasts |> gather_emmeans_draws()
mtsqrt2 %>% group_by(contrast) %>% dplyr::summarise(Prob = sum(.value>0)/n())

Summary plots

resp.plot <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=34.6, 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_fitted_draws(post.model1, 
                   n=500, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  ggplot(aes(x=Temperature, color=exp_group)) + 
  geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) + 
  geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) + 
  theme_classic(); resp.plot

Piecewise linear regression

Low-latitude - 28.5C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "28_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)

plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1            73   
## m = 2            64 83
## m = 3         55 70 85
## m = 4      40 55 70 85
## m = 5   15 40 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                 0.73     
## m = 2                 0.64 0.83
## m = 3            0.55 0.7  0.85
## m = 4        0.4 0.55 0.7  0.85
## m = 5   0.15 0.4 0.55 0.7  0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS    2.7785    0.5048    0.2062    0.1216    0.1174    0.1116
## BIC  -65.3272 -226.6618 -306.9730 -350.5902 -344.8844 -340.7792
plot(bp.data) 

breakpoint.values=c(54,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.86963 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.86963 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.86963 32.35335 0.03176730
## psi2.x 35.04506 33.76564 0.02974222
## psi3.x 36.14703 35.09270 0.02693903
## $x
##             Est.    St.Err. t value CI(95%).l CI(95%).u
## slope1 -0.019400 0.00076011 -25.522 -0.020909 -0.017890
## slope2  0.036286 0.00174150  20.836  0.032828  0.039745
## slope3  0.103900 0.00188890  55.003  0.100140  0.107650
## slope4  0.171910 0.00087702 196.010  0.170160  0.173650
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Low-latitude - 30C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "30_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3   18       70 85
## m = 4   15    55 70 85
## m = 5   15 30 55 70 85
## 
## Corresponding to breakdates:
##                               
## m = 1                     0.78
## m = 2                 0.7 0.85
## m = 3   0.18          0.7 0.85
## m = 4   0.15     0.55 0.7 0.85
## m = 5   0.15 0.3 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS    2.9454    0.5269    0.2322    0.1865    0.1643    0.1634
## BIC  -59.4938 -222.3828 -295.1253 -307.8064 -311.2763 -302.6371
plot(bp.data) 

breakpoint.values=c(54,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.86963 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.86963 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.86963 34.42398 0.07889169
## psi2.x 35.04506 35.07879 0.12326095
## psi3.x 36.14703 35.83758 0.14302773
## $x
##            Est.    St.Err. t value CI(95%).l CI(95%).u
## slope1 0.017059 0.00099279 17.1830  0.015088  0.019031
## slope2 0.107630 0.01762500  6.1065  0.072621  0.142630
## slope3 0.182670 0.01503100 12.1530  0.152810  0.212520
## slope4 0.238360 0.00529410 45.0230  0.227840  0.248870
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic() 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

Low-latitude - 31C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "31_core")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1            75   
## m = 2            66 84
## m = 3         55 70 85
## m = 4      40 55 70 85
## m = 5   15 40 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                 0.75     
## m = 2                 0.66 0.84
## m = 3            0.55 0.7  0.85
## m = 4        0.4 0.55 0.7  0.85
## m = 5   0.15 0.4 0.55 0.7  0.85
## 
## Fit:
##                                                                      
## m   0          1          2          3          4          5         
## RSS    1.90040    0.34415    0.14201    0.08957    0.08432    0.08360
## BIC -103.31265 -264.97683 -344.28312 -381.16143 -377.99555 -369.63983
plot(bp.data) 

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.94309 33.32652 0.02649379
## psi2.x 35.04506 34.50521 0.02054399
## psi3.x 36.14703 35.56366 0.02483205
## $x
##             Est.    St.Err. t value CI(95%).l CI(95%).u
## slope1 0.0065759 0.00028252  23.276 0.0060148  0.007137
## slope2 0.0473010 0.00137960  34.286 0.0445610  0.050041
## slope3 0.1127400 0.00152030  74.157 0.1097200  0.115760
## slope4 0.1620300 0.00079966 202.630 0.1604400  0.163620
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 28C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "28_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               79
## m = 2            70 85
## m = 3   15       70 85
## m = 4   15    55 70 85
## m = 5   15 40 55 70 85
## 
## Corresponding to breakdates:
##                               
## m = 1                     0.79
## m = 2                 0.7 0.85
## m = 3   0.15          0.7 0.85
## m = 4   0.15     0.55 0.7 0.85
## m = 5   0.15 0.4 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS    3.5327    0.5780    0.2508    0.2388    0.2338    0.2296
## BIC  -41.3137 -213.1309 -287.3916 -283.1020 -276.0267 -268.5972
plot(bp.data) 

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.94309 34.42427 0.07892526
## psi2.x 35.04506 35.07879 0.12333188
## psi3.x 36.14703 35.83765 0.14306228
## $x
##             Est.   St.Err. t value  CI(95%).l CI(95%).u
## slope1 0.0034725 0.0012763  2.7208 0.00093769 0.0060072
## slope2 0.1197800 0.0226570  5.2867 0.07478400 0.1647800
## slope3 0.2162000 0.0193220 11.1890 0.17782000 0.2545700
## slope4 0.2877800 0.0068057 42.2850 0.27426000 0.3012900
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 30C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "30_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3         55 70 85
## m = 4   15    55 70 85
## m = 5   15 37 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                      0.78
## m = 2                  0.7 0.85
## m = 3             0.55 0.7 0.85
## m = 4   0.15      0.55 0.7 0.85
## m = 5   0.15 0.37 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS    3.8094    0.6279    0.2542    0.2257    0.2236    0.2225
## BIC  -33.7721 -204.8428 -286.0614 -288.7430 -280.4479 -271.7575
plot(bp.data) 

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.94309 34.13055 0.04832474
## psi2.x 35.04506 34.92601 0.05935603
## psi3.x 36.14703 35.75668 0.07378457
## $x
##            Est.    St.Err. t value CI(95%).l CI(95%).u
## slope1 0.002638 0.00079262  3.3282 0.0010638 0.0042122
## slope2 0.100910 0.00938690 10.7500 0.0822690 0.1195600
## slope3 0.205170 0.00938690 21.8570 0.1865300 0.2238100
## slope4 0.279220 0.00354790 78.6990 0.2721700 0.2862700
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

High-latitude - 31C

posterior_draws <- lat_resp_dat2 |> 
  group_by(region, dev.temp,dev.temp.region) |> 
  data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature), 
                          to=max(lat_resp_dat2$Temperature), 
                          length.out=100), 
            mass=0.00148299) |> 
  add_linpred_draws(model1.gamk4, 
                   n=300, 
                   re_formula=NA) |>  
  unite("exp_group", c("region","dev.temp"), remove=FALSE) |> 
  filter(dev.temp.region == "31_leading")  

spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4) 

ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) + 
  geom_point(color="black", alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y, 
                                        color="red"), 
              size=2) 
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

strucchange analysis

spline.data <- as.data.frame(x= spline_fit$x) |> 
  mutate(y = spline_fit$y)|> 
  rename(x = `spline_fit$x`)
plot(spline.data)

fs.data <- Fstats(spline.data$y~1) 
bp.data <- breakpoints(spline.data$y~1)
summary(bp.data) 
## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
## 
## Breakpoints at observation number:
##                       
## m = 1               78
## m = 2            70 85
## m = 3         55 70 85
## m = 4      32 55 70 85
## m = 5   15 32 55 70 85
## 
## Corresponding to breakdates:
##                                
## m = 1                      0.78
## m = 2                  0.7 0.85
## m = 3             0.55 0.7 0.85
## m = 4        0.32 0.55 0.7 0.85
## m = 5   0.15 0.32 0.55 0.7 0.85
## 
## Fit:
##                                                                
## m   0         1         2         3         4         5        
## RSS    2.1641    0.3569    0.1421    0.1233    0.1224    0.1223
## BIC  -90.3197 -261.3436 -344.1864 -349.1651 -340.7615 -331.5675
plot(bp.data) 

breakpoint.values=c(55,70,85); spline.data[breakpoint.values, c("x")] 
## [1] 33.94309 35.04506 36.14703

segement analysis

Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints

my.lm <- lm(data=spline.data, y ~ x) 

spline.data[breakpoint.values, c("x")]
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm, 
                    seg.Z = ~x, 
                    psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg) 

my.seg$psi ; slope(my.seg)
##         Initial     Est.     St.Err
## psi1.x 33.94309 33.83815 0.03068568
## psi2.x 35.04506 34.78840 0.02982365
## psi3.x 36.14703 35.72150 0.03962484
## $x
##              Est.    St.Err.  t value  CI(95%).l  CI(95%).u
## slope1 -0.0032692 0.00039822  -8.2096 -0.0040601 -0.0024783
## slope2  0.0621610 0.00328720  18.9100  0.0556320  0.0686900
## slope3  0.1480500 0.00328720  45.0370  0.1415200  0.1545800
## slope4  0.2054600 0.00159820 128.5600  0.2022800  0.2086300
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)

ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()

ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
  geom_point(alpha=0.05, size=1) + 
  geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red", 
              size=3) + ylab("rate.output")+
  geom_vline(xintercept = my.seg$psi[1,2]) + 
   geom_vline(xintercept = my.seg$psi[2,2])+ 
  geom_vline(xintercept = my.seg$psi[3,2])+ 
   geom_vline(xintercept = my.seg$psi[1,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[2,1], color="blue") + 
  geom_vline(xintercept = my.seg$psi[3,1], color="blue") + 
  geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
  theme_classic()
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'